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An exact and explicit ladder-tensor-network ansatz is presented for nonequilibrium steady state 
of an anisotropic Heisenberg XXZ spin- 1/2 chain which is driven far from equilibrium with a pair 
of Lindblad operators acting on the edges of the chain only. We show that the steady-state density 
operator of a finite system of size n is - apart from a normalization constant - a polynomial of 
degree 2n — 2 in the coupling constant. Efficient computation of physical observables is facilitated in 
terms of a transfer-operator reminiscent of a classical Markov process. In the isotropic case we find 
cosine spin profiles, 1/n 2 scaling of the spin current, and long-range correlations in the steady state. 
This is a fully nonperturbative extension of a recent result [Phys. Rev. Lett. 106, 217206 (2011)]. 
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Introduction.- The Heisenberg model [TJ of coupled 
quantum spins 1/2 is the oldest many body quantum 
model of strong interactions. In spite of being ex- 
tremely simple it exhibits (in particular its anisotropic 
version, the XXZ model) a rich variety of equilibrium 
and nonequilibrium physical behaviors. In nature it pro- 
vides an excellent description of the so called spin-chain 
materials [5], and it is believed to provide the key for 
understanding of various collective quantum phenomena 
in low dimensional strongly interacting systems, such as 
magnetic or superconducting transitions in two dimen- 
sions. Although equilibrium (thermodynamic) proper- 
ties of XXZ chain are well understood in terms of Bethe 
Ansatz (BA) [3] , as the model represents a paradigmatic 
example of quantum integrable systems, its nonequilib- 
rium properties at finite temperature are subject to lively 
debate [4]. 

Ground states of strongly correlated systems generi- 
cally satisfy area laws [5] for block entropy characteriz- 
ing bipartite quantum entanglement, so they can be effi- 
ciently described by the so called matrix product states 
(MPS) or more general tensor networks [6j. MPS of 
small rank can provide even exact description of ground 
states, say in valence bond solids exemplified by the fa- 
mous AKLT model [7J. In fact, even BA eigenfunctions 
can be written in terms of MPS [5]. On the other hand, 
using the approach of open quantum systems and Marko- 
vian master equations [9], nonequilibrium steady states 
(NESS) of large one-dimensional locally interacting and 
dissipationless quantum systems put between a pair of 
unequal macroscopic reservoirs [THl E] , can be described 
in terms of a fixed point, or 'ground state' in the Li- 
ouville space, for a Hermitian super-operator with non- 
Hermitian boundary terms |12j . Application of the den- 
sity matrix renormalization group (DMRG) for simula- 
tion of such problems showed that a sort of super-area law 
is generically valid, and the density operator of NESS can 
be well described by a matrix product operator (MPO) 
of low rank [13] . However, no far-from-eqilibrium ana- 
logues of AKLT model have been known so far, and the 



purpose of this Letter is to show an explicit construction 
of an exact MPO form of NESS for a boundary driven 
XXZ spin chain. More precisely, a matrix element of the 
many-body density operator is a contraction of a very 
appealing ladder tensor network (LTN). 

We have recently proposed a new method [14] to solve 
for a Liouvillian fixed point of the XXZ chain, perturba- 
tively in the system-bath coupling constant. This method 
which expresses NESS in the form of a MPO with near- 
diagonal infinite rank matrices - reminiscent of a classi- 
cal Markov process in the auxiliary space - suggests new 
ways of integrability of strongly nonequilibrium quan- 
tum lattice gasses and appears to be unrelated to BA. 
In this Letter we show that - quite nontrivially - a fully 
nonperturbative extension of this method exists (in the 
strong driving limit of maximal bias, fx = 1 in nota- 
tion of Ref. [H]), with the constituent matrices satis- 
fying the same cubic matrix algebra (essentially different 
from quadratic algebras characterizing exactly solvable 
classical probabilistic lattice gasses, the so-called exclu- 
sion processes [15 ), but with modified boundary rela- 
tions. From our exact analysis, we: (i) prove ballis- 
tic transport (size n independent spin current) in the 
easy-plane regime, (ii) derive coupling independent co- 
sine spin-profiles, 1/n 2 scaling of the spin current and 
long-range spin-spin correlations in the isotropic regime, 
and (iii) prove insulating behavior in the easy-axis regime 
with kink-shaped spin profiles and exponentially (in n) 
decaying currents. We note that the physics of near- 
equilibrium XXZ chain is essentially different. There one 
has perturbative and numerical evidence of spin diffusion 
[TBI fT7] in the easy-axis regime, and alternative super- 
diffusive anomalous scaling in the isotropic point |16j in- 
dicating very rich phenomenology of the model. 

Nonequilibrium steady state.- We consider the Marko- 
vian master equation in the Lindblad form [91 111] 

^ = -i[H,p(t)}+J2^ kP (t)L{ {L{L k ,p(t)} (1) 
k 

for an open XXZ spin 1/2 chain with the Hamiltonian 
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H 



2afa j+1 +2a j af +1 + Aa*j<j*j +1 and 
symmetric Lindblad driving of coupling strength e acting 
on the edges of the chain only L\ — \fea\ , L 2 = y/ea~ . 
We write Pauli operators on a tensor product space 



Tn = (C 



2\®rc 



as a° 



Id being a 



d-dimensional unit matrix, where cr* = \(<J* ± i<r y ) and 
<7 x ' y ' z are the standard Pauli matrices. 

NESS is a fixed point of the flow ([I]) p^ = lim t _ i . 00 p(t) 

- i[H, poo] + eV Poo = 0, (2) 

with the dissipator map 



The quantum magnetic transport model (Tl I can be de- 
rived [18] by using standard assumptions 9 , or alter- 
natively, from an exact microscopic protocol of repeated 
interactions (19[ [20] where the left /right boundary spins 
are repeatedly and frequently put into arbitrarily strong 
contact with fresh up/down polarized magnets |21j . We 
shall now construct an explicit form of p^ in terms of 
the LTN ansatz, or equivalently, in terms of a product of 
two MPOs, which is exact for any value of the coupling 
parameter e. In fact, our simple explicit form allows us 
to study analytic dependence of NESS on e. 
Theorem. The normalized fixed-point solution of Eq. 
^ reads p^ = (tiR)~ 1 R with 

R = S n Sl (4) 

and S n a non-Hermitian matrix product operator 

S n = Yl (0\A Sl A S2 - ■■ A Sn \0)o- Sl <E>a S2 ■■■<g>a s " (5) 
(si,...,s„)e{+,-,o}' 1 

where a = I2 and Aq,A± is a triple of near-diagonal 
matrix operators acting on an infinite-dimensional aux- 
iliary Hilbert space % spanned by an ortho-normal basis 
{|0),|1),|2},...}: 



A = 



A + = k|0)<l|+5>+|r)<r+l|, 



(6) 



r=l 



A_ = \l)(0\+J2ar\r + l)(r 



with matrix elements (writing A := arccos A £ R U il 

-A) 
lA ' 

csin((2fc-l)A) sin(2fcA) 



n . sm(rA) 

a" = cos (rA) + ie— — — — , 
2 sin A 



a 2k-i = c sm (2fcA) + IE 
a 2k = csm (2fcA) — IE 



2(cos((2fc-l)A) + T 2fc _i)sinA 
c(cos (2kX) + r 2 k) 



' :2; 2sinA 

sin((2fc-l)A) . cos((2fc-l)A) +T 2fe _i 



(7) 



1 2fe-l 



IE 



2csin A 



sin((2Jfe+l)A) sin(2fcA)sin((2fc+l)A) 

a 2k = 1£ " 



2c(cos (2feA) + T2k) sin A 



Constant c e C — {0} and signs r r € {±1} are arbitrary, 
i.e. all choices of c, r r give identical operator S n 

Proof. We start by showing the following useful identity 

[H, S n ] = -is{a z ® 5 n _i - S n -i ® cr z ). (8) 

It is important to observe that the ansatz ^ does not 
contain any a z operator, while [H, S n ] can contain only 
terms with a single a z . Eq. (JsJ) implies that all the terms 
of [H, S n ] where cr| appear in the bulk 1 < j < n should 
vanish, resulting in exactly the same argument as in |14j 
leading to the same eight 3-point algebraic conditions: 



Vp := 2a+pa 1 -{a 1 , p} + 2a n pa+ - {a+a n ,p}. (3) [A ,A±A q 



0, 



{A ,A|} = 2AA ± A A±, 



{A T ,A^} 



Ai|-2A ± A T A±, 



2A{A2,A ± }-4A A±Ao 
2A[A2,A±] = [A t ,A 2 ± ]. (9) 

However, the boundary conditions should be different as 
in the perturbative case [13] ■ Namely the remaining set 
of terms where <r| appears at j = 1 or j = n in [H, S n ] 
is reproduced exactly by the right-hand-side of ^ if the 
following additional algebraic conditions are satisfied 

(0|A_ = (0|A + (A_A + - iel) = (0\A + A 2 _ = 0, 



A+|0) = (A_A+ - iel)A_|0) = A^A_|0> = 0, 



|A = (0|, A o |0) = |0), 



■IA, A J 



it. 



(10) 



Note the simple extra terms with amplitude — ie in com- 
parison to Eqs. (12) of Ref. [TJ]. Indeed, in order to get, 
e.g. a term with a\ in [hi, S n ], s± in ^ has to be + and 
then the condition (0|A + (A_A + — iel) = ensures that 
exactly S n -i will be constructed on the sites (2, . . . , n). 

Verifyin g ([9] ) and dXOft , which imply Q, for the repre- 
sentation ( 6|7 l results in trigonometric identities. 

The rest of the proof is to show that ([8J implies 
or i[H, R] = eDR (a). Left-hand-side of (a) can be 
transformed to [iH,S n ]Sl + S n [iH, S n }^ = e{S n (cr z ® 



S n (S n 



g z ) + (<r z ® s n -M - (S n 



°n-l 

ct z )S^} (b). Equation (10) implies that the first left- 



most (right-most) nontrivial operator of every term of 
S n is er + (c~). Thus we write S n =: a eg) S n -i + 
a + (g) P n -i =: S n -i ® cr° + Qn-i ® c~, defining Q n -i 
and P n _i as operators over T n -i, so the expression (b) 
further equals e{2a z ® <S' rl _iS' I t l _ 1 - cr+ (g) Pn-i^l-i - 



Sn— lP n _i 2S! n _iiS' ?1 _j 



^n-iQ^-i CS5 cr + } (c). On the other hand, writing the 
dissipator ^ as a sum of two local terms T> = X>l ® 
l„-i + ln-i <8> i>R, we have for the right-hand-side of 
(a): (eV h ® i„_!)(5 n 5t) + (i n _j ® e © R )(5 n 5j) (d). 
The first term of (d) can further be written out as 
(eP L ® 1«- 1 ) [(cr° ® 5„_ 1 + cr+ ® P„_ 1 ) (<r° <g> Sl_ 1 + a- ® 
Pt_ x )] = e P L (a°) ® S n _iS£_i +e2? L (a-)® ^-iP^i + 

£P L ((7+) ® Pn-l^Ui + £P L (CT+C7-) ® P^P^. SinCC, 

X> L (cr°) = 2cr z , Pl^) = -o- ± , P L (cr + cr _ ) = 0, we ar- 
rive at exactly the first three terms of (c). In an analo- 
gous way the second term of (d) results in the last three 
terms of (c). QED 
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FIG. 1. LTN contracting to a NESS density matrix element 
Thin (thick) lines represent bond dimension 2 (d). 



Corollaries.- Let us now derive some implications of 
our ansatz |4|5|: (i) Let \v),v= [y\,v%^ . . . , v n ) e {0, l} n 
denote the canonical many-body basis of J- ni oAy) = 
(1 — 2vj)\v). Then the matrix elements of MPO (15]) can 
±1) 



be written out as (± 

{l/\Sn\v) = 



A , A / • • • A 



jo). (11) 



(ii) This many-body matrix is upper triangular, i.e. 
(u'\S n \u) = if scq(V) > seq(i/) (where seq(i/) :— 
E"=i3- 2 " _i ) following from (0|A_ = dTob hence 
Eq. Q is the Cholesky decomposition of the many-body 
density matrix. We also have unit diagonal (v\S n \v) = 1, 
following from (0|A = (0| implying that NESS 

is always of full rank, (iii) Inserting the identity 1 — 
S/x l/ L< )(A < l in to B, the matrix elements of density oper- 
ator are obtained via contraction of a LTN (Fig. [lj 

(v'\R\k)= E (0|A m -,;A^- 4 ---A Mn _ < |0) 

pe{o,i}" 

x (0| A M1 _^ A M2 _, 2 ■ • • A^JO) . (12) 

A s denote the complex-conjugate matrices, obtained 
from ^ by complex-conjugating the amplitudes (J7|, 
equivalent to flipping the sign of e, A s — A s |_ e for a 
suitably chosen c (say as in Ref. [Sj). (iv) As the ma- 
trices ^ represent a nearest-neighbor hopping process 
in the auxiliary space 'H, they can - for any fixed chain 
length n - be truncated to a finite d = 1 + \n/2\ dimen- 
sional Hilbert space Wd spanned by {|0), |1) . . . \d — 1)}, 
still making the expressions ( |5|11|12[ ) exact, (v) bmcc 
hopping amplitudes (6[7| are all linear functions of the 
coupling e, the un-normalized NESS density operator R 
is a polynomial in e of degree not larger than 2n. In fact 
the degree is 2n — 2 as easily checked by explicit compu- 



tation, (vi) LTN (12 1 can be understood as an MPO on 
a tensor product auxiliary space 71®%, namely 

n 

R= £ (0|®<0|B Sl B S2 ---B s JO)«)|0)J]a-f , (13) 

se{o,±,z}" j=i 

introducing effectively d 2 x d 2 dimensional matrices 

B s = (WV)" 1 £ <',,A^®A^. (14) 
vy ,^e{o,i} 



Computation of observables.- Eq. (13) is a starting 
point for computation of expectations of physical observ- 
ables (A) = trpooA = trRA/trR. The normalization 
constant is computed as tri? = 2 n (0| ® (0|B "|0) ® |0), 
and a general expectation of a Pauli operator product 



reads (]\ j=1 <Tj - ,, llBi , „ 2 



nj =1 



For observables which are only products of <r|, say mag- 
netization profile (o~j), spin-spin correlations (ct|ct^), etc., 
one can use the same trick as in j!4j to further sim- 
plify the calculations. Namely, Bo and B z leave the 
auxiliary subspace /C of diagonal vectors, spanned by 
{|r) ® \r),r = 0, 1, 2, .. .}, invariant, B , z /C C /C. As the 
initial (final) vector |0)<g)|0) is also a member of IC, we can 
reduce the domain of our operators to K, defining tridi- 
agonal transfer matrices (TMs), T :— Bo|k;,V := B z |yc, 
or explicitly - using identification \r) ® \r) — > \r) 

oo i _|_ ,2 1—1^ 

T = E(l fl "| 2 I'M + ^V)(r+l\ + ^-Ir+lXrl), 



I + 1 I — I 

v = ^^|aj_| P><r+1 |_l«|L| r+ i ><r |) ) 



(15) 



r=0 



where we supplement ([7j) by a[] :— 1, :— ie, := 1, 
so that the physical observables are computed in terms 
of d x d matrix products 

(a]) = (0\T?- 1 yri n ~ i \Q)/{0\T n \0), (16) 
(a)a%) = (0|T'- 1 VT*- J '- 1 VT n -*|0)/(0|T n |0), etc. 

Another class of interesting physical observables are the 
spin current Jj = i(aj ~crj +1 — aj crT^) , local energy hj, or 
similar, which can be all formulated in terms of expec- 
tations of a non-Hermitian one-sided hopping operator 
wj := a~J The product B + B_ also leaves the diag- 
onal space K, invariant, so we introduce another vertex 
operator W := lB + B_|^, or explicitly 



(|a+|>)<r+l| + |a-|>+l)<r|) 
+ (a°) 2 a+a-|r)(r|+ a +a-(a0 +1 ) 2 |r+l)(r+l|}, (17) 



in terms of which the hopping expectation reads as 
( Wj ) = (0|T J '- 1 WT"- J '- 1 |0)/(0|T n |0). 



(18) 



Eqs. ( 15|17|7 ) imply ImW = — |T so the spin current 
(Jj) = — 2ha(Wj) is independent of the position j, man- 
ifesting local conservation law of magnetization. 

Let us now discuss some explicit results, graphically 
summarized in Fig. [2] We note that formulae ( 16|18 ) give 
efficient computational prescription which yields any ob- 
servable of this type in 0(n 2 ) arithmetic operations. In 
order to ensure numerical stability and to avoid singu- 
larities we suggest to choose the signs Tfc in computa- 
tion of auxiliary hopping amplitudes Q as = 1 for 
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FIG. 2. (color online). Spin profiles (uj) at n = 100 (a), and 
spin currents (J) vs. size n (b), for A = 3/2 (dashed), A = 1 
(dotted/blue), A = 1/2 (full curves), all for three different 
couplings e = 1, 1/5, 1/25 using thick, medium, thin curves, 
respectively. Red full curves show closed-form asymptotic re- 
sults [see text]: (cr|) = cos7ri^-, (J) = 7r 2 e 
in the main panels (a,b), and (J) oc e 



n arcoshA 



2 for A = 1 
in (b)-inset. 



cosfcA > 0, and = —1 for cosfcA < 0. For certain 
values of parameters even closed form results can be ob- 
tained. Analogously to the perturbative case [2], TMs 
have effective finite rank to + 1, i.e. they close on H m+ i, 
for a dense set of anisotropies, A = nl/m, which densely 
cover the easy-plane regime |A| < 1. This happens be- 
cause then a+ = for odd to, or = for even to, 
and the auxiliary hopping process gets cut. For example, 
for A = 1/2 = cos7r/3, we calculate spin profiles and 
currents by iterating a reduced TM T' = T|^ 3 



T' 



e 2 /2 
(l + e 2 )/4 







)/24 



(9 

3(l + e 2 )/8 (l + e 2 )/4 



(19) 



V = V|„ 3 ,W 



combined with the reduced 3x3 vertex matrices 
W|ft 3 , Explicit expressions for 
fe /, \Jj) can easily be obtained by means of di- 
agonalization of T'. We obtain exponential convergence 
towards the thermodynamic limit (TL), n — ¥ oo, with the 
rate given by the ratio of two leading eigenvalues of T', 
and asymptotically flat spin profiles (a 7 -) ps (Fig. [2^,). 
We prove ballistic transport by explicitly computing the 

r V/T\l (V81+74 £ 2 +9£ 4_ 7 _ 3e 2- )e n 

limit (J,-)|„^oo = 4(1+e 2) — (Fig. pb), hav- 
ing a non-monotonic e-dependence starting as ~ e/2 for 
small e (consistent with [14), a maximum at e* ps 1.63, 
and decaying asymptotically as ~ 4/(3e) for large e, 
qualitatively agreeing with similar results for the non- 
interacting XX [19] and XY chains [22]. Similar finite 
dimension analysis can be made for some larger denomi- 
nators m. On the other hand, for A > 1, the TM T has 
always an infinite rank. In the easy-axis regime |A| > 1, 
explicit computations reveal almost e-independent kink- 
shaped spin density profile (Fig. [2^) - agreeing with nu- 
merical simulations of negative differential conductance 
[IB] - and asymptotically exponentially decaying current, 
(Jj) oc (IAI+VA 2 - I)"" (Fig.^-inset), consistent with 
suggested ideally insulating behavior [2"3"] . 



At the end, let us briefly focus on the isotropic case 
A = 1. In this case, our hopping matrices Q have to 
be regularized by taking T2k-i = 1, T2fe = — 1, and c = 
1/A before taking the limit A — > 0, yielding the hopping 



amplitudes: a° = l+kr/2, Oq = is, 



+ 

2/c-l 



= 2fc+iefc(fc- 



= 2k+iek 2 



«n 



l)l a 2k ~ " 1TKIV ' "0 — x > "2fc-l — lc > "2fe — 1 2 

The following formulae can be verified with some effort 



= 1, a,,, , = ie, a ,„ = ie(k+\)/k. 



[T, [T,V]] 

(0|(T-V) 
(0|T n |0) 



£ 



(2V + {T,V}) ; 

(T + V)|0) = 

(4n - 3) 2 



10), 



(20) 
(21) 



32tt 2 



-a +l + 0{n- L ), (22) 



where a w 0.0346. Multiplying (J20j) by (OlT-?" 1 from 
the left, and T""^ 2 |0) from the right, and using {22} 
we obtain in the continuum limit M(x = j^Ej) '■= (o~ 
a differential equation M"(x) = -tt 2 M(x) + O(^) 



3 ' 

and 



from (21) the boundary conditions M(0) = — M(l) = 
0(— ), yielding a magnetization profile 



cos7ra;, or <cr" 



cos7T^3Y, for arbitrary e ^> 



1 

M(x) 

e* = 27r/n (Fig. [2^). "Similarly we use ( |20p2| and the 

continuum approximation to calculate the connected cor- 
= ii-i „, — fc- 



relator C(a; 



,2/ 



3) 



^)-(apK), for 



j 7^ fc. However, as it turns out that the leading or- 
der O(n ) of C(x,y) exactly vanishes, we solve the cor- 
responding differential equations perturbatively in the 
next order in 1/n. Straightforward but tedious calcula- 
tion gives C(x,y) ~ ^/(min(ar, y), m&x(x, y)) + 0{^), 
where f(x, y) = 2irx(y—l) sin(7r:r) sin(7ry) +cos(7rx)((l — 
2y) sin(7ry) +ir(y~ l)y cos(7ry)). This is another, now an- 
alytic, indication of long-range correlations in far from 
equilibrium quantum NESS recently observed numeri- 
cally or in non- interacting systems [24]. Eq. (22) and 
ImW = — §T imply anomalous sub-diffusive scaling 
(Jj) = §(0|T n - 1 |0)/(0|T n |0) fa Tr^n- 2 , again valid 
for any e 3> e* (n) (Fig.J^)). For s -C e* we reproduce the 
perturbative result [H], (Jj) = \e, (of) = \e 2 (rt+l-2j) . 

Discussion.- An explicit LTN/MPO ansatz has been 
written describing the many-body density matrix of 
NESS of strongly boundary driven XXZ chain, for any 
bath-coupling strength. Computation of the physical ob- 
servables in NESS is facilitated in terms of tridiagonal 
transfer matrices which are reminiscent - except for non- 
conservation of 'probability' - of a classical Markov pro- 
cess in the auxiliary space. Results in TL can be obtained 
by studying the spectral properties of the transfer opera- 
tor. Studying Liouvillian gap or relaxation rates to NESS 
and related uniqueness of NESS is yet to be addressed. 
It also remains open to what extend our solution (|4][7} 
can be generalized to other bath-models, for example the 
case of weak driving has fundamentally different physi- 
cal properties [13l [16] . Our method seems to open a new 
ground for constructing exactly solvable nonequilibrium 
quantum problems in one dimension, and seems to be 
unrelated [35] to existing algebraic methods [HI [T5] ■ New 
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exactly solvable models could perhaps be constructed by 
studying alternative cubic algebras of type ^ . 

Discussions with M. Znidaric and support by the 
grants J 1-2208 and P 1-0044 of ARRS (Slovenia) are ac- 
knowledged. 
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